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Abstract: We construct an approximate expression for the total cross section for the 
production of a heavy quark-antiquark pair in hadronic collisions at next-to-next-to-next- 
to-leading order (N 3 LO) in ot s . We use a technique which exploits the analyticity of the 
Mellin space cross section, and the information on its singularity structure coming from 
large N (soft gluon, Sudakov) and small N (high energy, BFKL) all order resummations, 
previously introduced and used in the case of Higgs production. We validate our method by 
comparing to available exact results up to NNLO. We find that N 3 LO corrections increase 
the predicted top pair cross section at the LHC by about 4% over the NNLO. 
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1 Introduction 

The lack of discovery of any new physics signal so far has made of the Large Hadron Collider 
(LHC) even more a precision machine than it ever was. Because the LHC is a hadron 
collider, the largest uncertainties are related to nucleon structure, i.e. parton distributions, 
and higher order QCD corrections. This has consequently led to enormous progress in the 
computation of higher order corrections to QCD processes in the last several years, recently 
leading to the publication of first N 3 LO results for a hadron collider process [1]. 

In Ref. [2] some of us have proposed a general methodology for the determination of 
approximate expressions for higher-order corrections to QCD processes. The basic idea is 
that QCD corrections to many important hard processes are known to all orders in the 
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strong coupling in two opposite limits: the high energy limit, in which the available center- 
of-mass energy is much greater than the invariant mass of the final state, and the soft limit, 
in which the invariant mass is close to threshold. By performing a Mellin transform, this 
knowledge can be turned into information on the singularities of the hard (partonic) Mellin- 
space cross section to any given finite order, viewed as an analytic function of the Mellin 
complex variable: the soft limit determines the singularity at infinity, and the high-energy 
limit determines the rightmost pole on the real axis. It is then possible to reconstruct an 
approximate form of the function by exploiting this knowledge: indeed, knowledge of all 
the singularities would determine the function completely. 

The case of Higgs production, to which the methodology of Ref. [2] was first applied [2- 
4], is particularly simple in many respects: the leading-order process has fixed kinematics 
(only one scalar colorless particle in the final state), and the total cross section is known [5] 
to be dominated by its threshold limit. The validation of the methodology which comes 
from checking that it provides a good approximation to known results in the case of Higgs 
is thus not necessarily the most compelling. 

In this paper we turn our attention to heavy quark production: a process which is rather 
more complicated in terms of kinematics and color structure, with the goal of performing 
a more stringent test of our methodology. Next-to-leading (NLO) QCD corrections to this 
process were computed long ago [6, 7], though fully analytic expressions only became avail¬ 
able much more recently [8], and were a first necessary step towards the full determination 
of the NNLO result which was very recently achieved [9-12], These results revealed that 
the physical-space partonic cross section has a particularly complex singularity structure, 
which raises the question of how this relates to the Mellin-space singularities which are 
used for the approximation of Ref. [2], and also, that existing attempts at an approximate 
NNLO determination [13] were rather off the mark. 

Here, we will present an approximate determination of the N 3 LO QCD corrections 
to the total cross section for heavy quark production. In Section 2 we will review the 
singularity structure of the heavy quark production partonic cross section both in physical 
and in Mellin space. Specifically, we will discuss the relationship between the total cross 
section and the invariant mass distribution, to which the methods of Ref. [2] apply, and show 
that the peculiar singularity structure of the physical-space cross section leads nevertheless 
to the standard Mellin-space singularities. The rest of our treatment closely follows that 
of Higgs production: in Sects. 3 and 4 respectively we present the derivation of the soft 
and high-energy singularities from the respective resunnnations. They are combined in 
Section 5 where we present our final result: first, we compare the approximation obtained 
through our procedure to the known exact result up to NNLO, and then we present our 
final approximate N 3 LO both at the parton and the hadron level. 

2 Analytic structure of the partonic cross section 

The analytic structure of the cross section for the production of a heavy quark pair poten¬ 
tially differs from that for Higgs production discussed in Ref. [2] for two different reasons, 
thereby potentially hampering, or requiring some adaptation of the methods used in that 
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reference. First, the known [8] fact that the partonic cross section for heavy quark pro¬ 
duction already at NLO has unphysical singularities in the complex plane of its kinematic 
variables may suggest that the singularity structure of its Mellin transform also does not 
follow the general pattern discussed in Ref. [2], Second, as we shall see in more detail below, 
the analytic structure discussed in Ref. [2] is generic for partonic differential cross sections, 
while the case we are presently considering is that of a total cross section. We will see 
that the first issue is of no concern, while the second requires a careful discussion of the 
relation of the total heavy quark production cross section to the corresponding invariant 
mass distribution. 

We will address the two issues in turn, after a quick review of standard notation and 
general results. 


2.1 Notation and general results 


The total cross section for the production of a heavy quark-antiquark pair can be written 
in factorized form as 


v(m 2 ,p h ) = Ph^2 [ ~^ij 
ij P 




1 ~ ( 2 / 2 \ m m \ 

-Oij [m ,p, a s (pj, -w, -o- 
P \ Pf pi 


( 2 . 1 ) 


where ph = Am 2 / s , s is the hadronic center of mass energy, J£ij(p,p 2 ) are parton lumi¬ 
nosities, m is the mass of the heavy quark, /x F ,/% are factorization and renormalization 
scales, and the scaling variable is p = 4 m 2 /s, with s the partonic center of mass energy. 
To simplify notations, in the following we will not indicate explicitly the dependence of 
parton luminosities on the factorization scale p, F and that of the partonic cross section on 
the strong coupling and on p, F ,p R . 

The partonic cross sections aij admit a perturbative expansion in powers of the QCD 
coupling: 


&ij{m 2 ,p) = a^\p) + a s &lf(p) + a 2 a[f (p) + a 3 s &lf (p) + 0(a 4 s ) 


( 2 . 2 ) 


where, thanks to the m 2 prefactor, the coefficients &ij\p) are dimensionless. Eq. (2.1) is 
in the form of a convolution product, so its Mellin transform 

a(m 2 , N) = [ dp h p%~ 2 a(m 2 ,p h ) (2.3) 

Jo 

factorizes in terms of the Mellin space luminosity and partonic cross section function, defined 
respectively as 


&(N)= C dpp N - l J?(p) 

Jo 

(2.4) 

dij(m 2 ,N)= / dp p N ~ 2 dij(m 2 , p), 

Jo 

(2.5) 

a(m 2 ,N) = '^2^ , ij(N)aij(m 2 ,N). 

(2.6) 


ij 
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If the Mellin transform integral has a finite convergence abscissa, the JV-space partonic 
cross section is an analytic function of the complex variable N. dehned by the integral rep¬ 
resentation Eq. (2.5) to the right of the convergence abscissa, and by analytic continuation 
elsewhere. Therefore, it is fully determined by the knowledge of its singularities. 

In this paper, we will concentrate on the gg partonic channel, which is the most relevant 
(in the MS scheme) at LHC energies, while we leave the study of other parton subprocesses 
for future work. For this reason in the following we will drop the parton indices ij, with 
the understanding that i = j = g in all parton luminosities and partonic cross sections. 

The singularity structure of a generic differential partonic cross section is relatively 
simple. The singularity at infinity is determined by soft gluon radiation, which leads to 
a growth of the cross section with increasingly high powers of In N at higher perturbative 
orders. For finite N, singularities away from the real axis are not allowed, as they would 
have to come in pairs and would thus lead to oscillatory behaviour of the cross section at 
high energy. From Regge theory one expects that the rightmost singularity on the real axis 
is a multiple pole located at N = 1 [14], with further multiple poles along the real axis at 
N = 0,-1, —2,.... This expectation is confirmed by explicit fixed-order calculations. 

While knowledge of the residues of all poles is required in order to fully determine the 
function <r(m 2 , N), its behaviour in the region 1 < Re N < oo is mostly controlled by the 
rightmost pole at N = 1, with poles further to the left having an increasingly small impact. 
Using the saddle point approximation one can show that the hadronic cross section is mostly 
determined by the behaviour of the partonic cross section &(m 2 , N ) in the vicinity of a single 
(saddle) value of N on the real axis, to the right of the rightmost singularity [2, 15] (see 
also Sect. 5.1). In Ref. [2] it was suggested that knowledge of the rightmost singularity on 
the real axis and of the singularity at infinity is sufficient to determine the partonic cross 
section in the region where the saddle-point is located with reasonable accuracy. 

As mentioned in the introduction, however, in the specific case of heavy quark pro¬ 
duction, we encounter two difficulties. The first, is related to the fact that the physical 
momentum-space coefficient function has unphysical singularities. The second is due to the 
fact that the partonic cross sections &ij{m 2 , N) Eq. (2.6) vanish as N —> oo, rather than 
growing logarithmically. 

We will address both issues in turn. The first issue turns out to be of no concern: 
the analytic structure in Mellin space is unaffected by the unphysical momentum-space 
singularities. The second issue instead is due to the fact that the reason why the partonic 
cross sections discussed in Ref. [2] grow as N —> oo is that they are distributions, rather 
than ordinary functions: indeed, elementary properties of Laplace transforms imply that 
the Mellin transform of an ordinary function, if it exists, must vanish as N —> oo. Now, the 
integral of a distribution is an ordinary function, so it is clear that the behaviour of Ref. [2] 
can only hold for partonic cross sections which are sufficiently differential. As we shall see 
below, it is the invariant-mass distribution which behaves at large N in the way discussed 
in Ref. [2]. However, by exploiting the relation between total cross section and invariant- 
mass distribution, it is possible to relate their respective large-iV behaviors, and define a 
coefficient function whose singularity structure follows the general pattern discussed above. 
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2.2 Unphysical singularities in momentum space 

The perturbative coefficients (p) display a class of spurious singularities in the complex 
p plane outside the physical region 0 < p < 1. However, after Mellin transformation, even 
in the presence of such spurious singularities in p space, the result has the expected analytic 
structure, as we now show. 

As an example, we present the case of NLO corrections of top pair production in the gg 
channel. In Ref. [8], the complete analytic form of dgg (p) is computed, and it contains, in 
addition to the usual (physical) singularities in p = 1 and p = 0, four further singularities 
(branch points) located at p = —4, p = 4, p = —1, p = —^ (see in particular the functions 
Fi(p), i = 1,..., 4 of section 4 and the discussion of section 5 of Ref. [8]). For simplicity 
we focus our attention only on the first spurious singularity; similar considerations hold for 
all the others. 

The function 


f 2 (p) 


x{p) 


>x{p) 


dzf 2 {z ) 


(2 z + 3) ^12 In z In + 3 In 2 z + 12Li2 (—z) + i r 2 j 


4 (z 2 + 3z + 1) 


i -yr^-p 

1 fv'W 


(2.7) 

( 2 . 8 ) 

(2.9) 


has a logarithmic branch cut starting at p = —4, because the integrand has a simple pole 
in 

z = 3 = x(-4). (2.10) 

However, the Mellin transform is 


F 2 (N) 



( 2 . 11 ) 


the p integral is convergent for Re N > 1, and the result 

aN-i r i 

f 2 (AO = j- J o dzz N -\l + zf- 2N h(z) (2.12) 

can be analytically continued to the whole complex N plane, except the isolated points 
N = 1,0, —1, —2,..., where the result has simple poles. This can be seen for example by 
repeatedly integrating by parts using z^” 1 as a differential factor. 

This is the analytic structure one expects for the Mellin transform of a physical cross 
section. We conclude that the presences of unphysical singularities in the p space does not 
affect the general structure the partonic cross section in N space. Of course, it could be 
that the presence of new structures in the partonic cross section at higher perturbative 
orders affects the numerical size of the residues of singularities, for the same reasons why 
it makes standard scale-variation estimates of higher order terms unreliable. 
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2.3 Total cross section and invariant-mass distribution 

We will first show that the invariant-mass distribution for heavy-quark production has the 
large N singularity structure discussed in Ref. [2] and dominated by Sudakov radiation, 
then we will prove that the coefficient function C(N), defined by factoring out the leading 
order contribution in the total cross section: 


a(m 2 , N) = a LO (m 2 , N)C(N), 


(2.13) 


exhibits the same Sudakov enhancement at large N. 

The invariant-mass distribution £(m 2 ,£,.s) is defined as 

^{m 2 ,C,z) = s-^^(m 2 ,£,z), (2.14) 

which is a function of two dimensionless ratios, due to the presence of an extra energy scale 
M 2 , the invariant mass of the quark-antiquark pair. We choose to express it as a function 
of z = M 2 /s, and of £ = 4 m 2 /M 2 = p/z, and we insert the factor of s in Eq. (2.14) for 
later convenience. The total cross section is related to the invariant-mass distribution by 



(2.15) 


In Mellin space, the relation between a and £ can be obtained by computing the Mellin 
transform of Eq. (2.15): 


»i r i 

N -2 / 


a(m 2 ,N)= / dp p 

Jo 

= J dz J dpp N ~ 2 Tj (jn 2 , z'j . 


z / 

(2.16) 

z / 

(2.17) 


'o J o 

Changing integration variable we obtain 

d(m 2 ,N)= [ dzz N ~ 1 [ dZZ N - 2 X(m 2 ,£,z) = £(m 2 ,JV- 1,N). 
Jo Jo 


(2.18) 


Using this last relation, we will show that C(N ) implicitly defined in Eq. (2.13) has 
the same Sudakov singularities in the large N limit as the invariant mass distribution 
£(m 2 ,^, N) in the limit £ —» 1. This limit is called the absolute threshold limit, and it 
corresponds to the double limit £ — > 1 and z — > 1 of £ (m 2 , £, z) in z space. Indeed, we see 
from Eq. (2.15) that, when p —» 1, only the region z —> 1 (and hence £ —» 1) contributes to 
the integral. 

The behaviour of £(m 2 ,£, 2 :) for z close to one is governed by Sudakov resummation. 
In analogy with Drell-Yan or resonant Higgs production, the G(a™) perturbative coefficient 
for £(m 2 ,£, z) is a linear combination of the distributions 


TJ>k{z) 


ln fc (l - z) 


6 - 


1 - z 


<5(1 - z) 


(2.19) 







with 0 < k < 2n— 1, plus contributions which are less singular in the threshold limit z —>■ 1. 
Correspondingly, the perturbative coefficients of 

E(m 2 , £, N) = f dzz N ~ 1 i:(m 2 ,Z,z) (2.20) 

Jo 

grow at large N as ln fc N, 0 < k < 2 n. 

The resummation of these large logarithmic contributions has been performed in re¬ 
cent years, both in momentum space with Soft Collinear Effective Theory (SCET) tech¬ 
niques [16-18] and in Mellin space [19-21], In the latter case one finds [20] 


£ res (m 2 , £, N) = exp [G (N)] Tr i H(£, a s {M 2 )) 


X exp 


x exp 


IJo 


»i ~N—1 _ i 

dz —r s t (C,o;< j (M 2 (l - z) 2 )) 

“1 ~N—1 1 

dz Z ~ 


S [ a, [ 


1-z 


-T s (£,a s (M 2 (l-z) 2 )) 


( 2 . 21 ) 


up to corrections which vanish in the large -N limit. Here, bold symbols are used for matrices 
in color space. In the case of heavy quark pair production in gluon-gluon fusion, these are 
3x3 matrices because there are three independent colour configurations [21]. The function 
G(N ) is given by 


G(N) = 



z 


N -1 


- 1 


1 - 2 


r M 2 (l-zf j 2 

2/ ~^A{a s ^ 2 ))+D{a s {M 2 { 1 

hi h 



( 2 . 22 ) 


where A(a s ) [22, 23] and D(a s ) [24] have perturbative expansions in powers of a s , which 
are known up to O [oig) (see for instance [25-27]). The soft anomalous dimension Ts(£, a s ) 
and the soft (matrix) function S(a s ) originate from soft emissions in the presence of a 
heavy quark pair. Finally, the matrix function H(£, a s ) represents the hard contribution to 
the cross section. NNLL accuracy is achieved by expanding the cusp anomalous dimension 
A up to three loops, the characteristic function D and soft anomalous dimension Ts up 
to two loops, and the soft and hard functions S,H to one loop. Moreover to predict the 
term proportional to the delta function 5(1 — z) at 0(a 2 ), we just need the sum of the two 
loop contributions to the hard and soft functions, H and S, which can be determined by 
matching with NNLO calculation. 

However, what we are interested in is the absolute threshold limit £ —> 1 of Eq. (2.21). 
We now show that in this limit Eq. (2.21) greatly simplifies. For £ ~ 1 we can replace 
M 2 by 4m 2 in the argument of a s , the difference being suppressed by powers of 1 — £. 
Furthermore, it was noted in Ref. [28] that in this limit the matrix r§(£,a s ) is diagonal in 
the singlet-octet basis defined e.g. in Ref. [21]. To order a 2 it takes the form [28] 


Ts(£, Oi s ) = —C A 


a s / a s \ 2 , 

2WU (*+&-!) 


ns = r s (a s )n 8 , 


(2.23) 
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where 


K = I — - —) C* - —n,, 
'18 6 / 18 * 


(2.24) 


and the matrix II 8 is the projector over the octet subspace: 

/o 0 o' 


n 


8 — 


0 1 0 
0 0 1 


(2.25) 


We now turn to the soft matrix S(a s ). By a suitable definition of H, it can be chosen 
to be the unity matrix I at leading order. Its NLO expansion, sufficient to achieve NNLL 
accuracy, is given by 


S(a s ) = 1 + — a,II 8 + 0(a 2 ). 

7T 

Thus, in the absolute threshold limit, Eq. (2.21) reduces to 


(2.26) 


£ res (m , £, N) = exp [G(N)] Tr { H(£, a s (m 2 )) 


C a f 4m , „ 

1 + (^2“ ) n 8 + C>(a 


x exp 


,i n -l _ x 

2n 8 / dz— - 

Jo 1 - z 


Ts(a s (4m 2 (l - z) 2 )) 


(2.27) 


The calculation of the color trace in Eq. (2.27) is greatly simplified by noting that, for any 
complex number a, 

er - = e iig + iii 

where 


0 an s = e a ng + ni (2.28) 


iii + ii8 = i; iii 2 = Hi; n 8 2 = n 8 ; nin 8 = o. (2.29) 

As a consequence, the resummed invariant mass distribution splits into the sum of an octet 
and a singlet component: 


£ res (m 2 ,£, N) = exp [G(N)] Tr [H(£, a,(m 2 ))n g 


C A ( 4 m 2 \ | ^ 

1 H - ol & 

TV 


N 2 


+ 0{a 2 s 


x exp 


' f 1 , z N ~ 1 - 1 
2 dz 


T s (a s (4m 2 (l - z) 2 )) 


1 0 1 -z 

+ exp [G(N)j Tr [H(£, a,(m 2 ))n 1 ] [l + 0(a 2 )] . 


(2.30) 


The matrix H(£, a s ) cannot be expanded in powers of £ — 1, because it is proportional 
to the phase-space factor \Jl — £. However, one may factorize the leading order coefficient 
as in Ref. [28], to obtain 


H(£, a s ) = H(°)(0 [l + H«(£)«, + lK 2 )(£)a 2 + 0(a 3 s ) 
The coefficients H^(£), i > 1 are now analytic around £ = 1. Thus 

H(£, a s ) = H(°>(£) [l + h«a s + h^a 2 + 0(a 3 s ) + 0(1 - £) 


(2.31) 


(2.32) 




















where = H^'(l) are constant matrices. A further simplification arises from the partic¬ 
ular structure of the matrix Hi°'(£, a s ). Indeed, one finds 


H (0) (0« s ) = 


( o\°(m 2 ,£) a\°(m 2 ,^) 

2 a 5 -lo 


3 a\°(m 2 ,0 I 

V 0 


0 

2 o 

0 a£°(m 2 ,€) - \a\°{m 2 ,i) 


(2.33) 


where a\°{m 2 , £) are the leading-order total cross sections in each color configuration. Since 


g|°(m 2 ,Q _ 5 
?™<7L°(m 2 ,0 2’ 


G 


LO/ 


(m ,£) = o-L°( m ,£) + ££ ( m ,£) 


we have 


H(°)(£, a 8 ) = d LO (m 2 , Oh (0) + 0(1 - £); 


It follows that 


h(°) 


/? ?o\ 
f f 0 
y o o oy 


Tr [H(^,a s )ni] =«r LO (m 2 ,05i(a s ) +0(1-0; 1 = 1,8, 


(2.34) 


(2.35) 


(2.36) 


where 


gi(a s ) = Tr 


h (0) (l + a s h^ + a 2 h (2 ^ + 0(a®)) ni 



(2.37) 


are independent of N. 

Eq. (2.30) may be cast in a more familiar form, by reabsorbing the logarithmic terms in 
the factor l + C' A Q!<j(4m 2 /A r2 )/7r in a redefinition of the function This also generates 

IV-independent terms, which can be absorbed in a redefinition of the constants gg' 1 [28]. 
Putting everything together, we obtain 


E res (m 2 , £, N) = <r LO (m 2 , £) gi(a s ) exp [Gi(N)] + O (2.38) 

1=1.8 ' ' 


with 


r 1 Z N ~ 1 - i 
Gi(N) = dz 
Jo 


r4m 2 (l-z) 2 j 2 

2/ —y- A(a s (^ 2 )) + T> I (a s (4m 2 (l - z) 2 )) (2.39) 


'Mp 


1 - z 
T>i(a s ) = D{a s ) 

D 8 (a s ) = D(a s ) + 2Tg(a s ) + 2C A /3 0 (g)' + 0(a 3 s ). 


(2.40) 

(2.41) 


Finally, recalling the relation Eq. (2.18) between total cross section and invariant-mass 
distribution, we get 


d(m 2 ,N) = I dU N ~ 2 ^(rn 2 ,Z,N) 
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= a h °(m 2 , N) Y 9i{a s ) exp [Gi{N)} + O f (2.42) 

1 = 1,8 ' ' 

where 

<r LO (m 2 , N) = f d^ N - 2 d h °(m 2 ,t), (2.43) 

Jo 

is the Mellin transform of the Born level total cross section. By comparing Eq. (2.42) with 
the definition of the coefficient function Eq. (2.13), we see that 

C(N) = Y 9 i(a s ) exp [G^N)} + O Q-) (2.44) 

1 = 1,8 ' ' 

has the same singularity structure as the Mellin-transformed invariant-mass distribution in 
the N —>• oo limit. 

2.4 Analytic structure of the coefficient function 

We now focus on the coefficient function C(N), implicitly defined in Eq. (2.13). We have 
shown in the previous section that in the N —> oo limit C(N ) has the same singularity 
structure as the Mellin-transformed invariant mass distribution, whose behaviour in turn 
is determined by soft-gluon emissions. On the other hand, we note that C ( N ) has the 
same leading singularity in IV = 1 as a(m 2 ,N), because <r LO (m 2 ,lV) is subleading in the 
high-energy regime. Therefore, we can construct an approximation to C(N) according to 
the procedure of Ref. [2], as 


C appr ox(A0 = C soh (N) + C'h.e.(iV), (2.45) 

where C so ft contains the terms predicted by Sudakov (soft) resummation, and C'h.e. the terms 
predicted by BFKL (high-energy) resummation. Explicit expressions for both components 
are given in the following Sections. 

3 Large-iV contributions 

3.1 Threshold resummation and analyticity 

In this subsection, we will extend the procedure outlined in Ref. [2] for the Higgs production 
cross section to the coefficient function C(N) of Eq. (2.44). The cusp anomalous dimension 
A(a s ) has been computed up to three loops [23], and the colored characteristic anomalous 
dimensions Di(a s ) are known completely at two loops [29], allowing us to achieve NNLL 
accuracy. Hence, we are able to determine all terms of order a™ In'" N, with 2n — 3 < 
m < 2n, in the coefficients of the perturbative expansions of C(N). The inclusion of a 2 
contribution in the constant terms g\ enables the extension of our prediction to 2n — 4 < 
m < 2 n. We are thus able to predict all large-A non-vanishing contributions to C(N ) 
up to O (as), and all logarithmically enhanced contributions except the single log and the 
constants at O(a^). The single log at 0(af), formally a N 3 LL contribution, is produced 
by the third order of the colored characteristic anomalous dimensions D\. At this order 
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only the singlet is fully known [27], but we do not know the contributions to D g ' coming 

from the N'^LO soft anomalous dimension and soft function. We thus set and 

-13) 8 

< 7 j = 0 in the following. The extra uncertainty induced by this assumption on our results 

will be discussed in Sect. 5.3 below. 

As pointed out in Ref. [2, 5], the quality of the soft approximation to the full cross 
section significantly depends on the choice of subleading terms which are included in the 
resummed result. Since the resummation procedure only fixes the coefficients of logarith¬ 
mically divergent terms, there is a certain freedom in defining how the soft approximation 
is constructed, by including contributions that are suppressed in the limit N —> oo. The 
idea proposed in Ref. [2, 3] is to include subleading terms that are known to be present in 
the exact calculations, in order to preserve the known analytic structure at small N. 

In principle, the exponent Gi(N) Eq. (2.39) is ill-defined, because the integration range 
includes the Landau pole of the strong coupling. This problem is usually avoided by com¬ 
puting the integral in the large -N limit: the resummed IV-space result is then well-defined 
as a function of In N. The logarithm, however, has an unphysical cut at N = 0, where the 
cross section should only have poles. 

However, as pointed out in Ref. [2], if we are only interested in the expansion of 
Eq. (2.39) in powers of a s (m 2 ) to finite order, the problem of the Landau pole does not 
arise. The Mellin transform in Eq. (2.39) may then be computed exactly, provided the 
integrand is understood to be expanded in powers of a s (m 2 ) to some finite order. The 
result is a function of N which has the correct logarithmic behaviour at N —> oo, but is 
free of unphysical cuts on the real negative axis from N = 0. 

An explicit calculation yields 


Gi(JV) = J>? 

77.—1 


n 7i—l 


_k =o 


k =0 


where T>k(N) are the Mellin transforms of the distributions 


T>k{z) 


ln fc (l -z) 
1-z 


(3.1) 


(3.2) 


and the coefficients up to order n = 3 are given in the Appendix A, Eqs. (A.l) 

and (A.2), and depend respectively on A(a s ) and Di(a s ) only. Explicit forms of the func¬ 
tions T>^{N) are given e.g. in Ref. [2]; they are seen to only have poles on the real axis. 

Following the argument of Ref. [2], we now observe that the logarithmic enhancement at 
threshold arises from the integration over the transverse momentum of the emitted gluons, 
which has the form 


rM 


(!--) 


V 99 {z) 


dk'Y 


A g(z) 
1 - 2 


In 


1 - 2 


+ In 


M 

T 


(3.3) 


where V gg (z) is the gluon-gluon splitting function for z < 1 and A g {z) = (1 — z)V gg (z), with 
4,(1) = Ad) (see Eq. (A. 3)). Thus, it is natural to include subleading terms in Eq. (2.39) 
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by restoring the factor of 1 /y/z in the upper integration bound, and by replacing A ^ by 
the expansion of A g (z) about z = 1 up to some finite order (keeping the full expression 
of A g (z ) is not advisable, because an unphysical singularity in z = 0 would appear; see 
Ref. [2]). 

We therefore replace Gi(N) with 


Gi(N) = dz 


+ ci, 


v JV-l 


- 1 


1 - z 


'Mf 


d/i 2 


zA{a s {iJ 2 )) + Di ( a 



which differs from Gj(N) by non-logarithmically enhanced terms. The factor of z in the 
first terms arises after expansion of A g (z) in powers of 1 — z to first order; indeed, it turns 
out that 

Ag(z) = Ag( 1) - (1 - z)Ag{ 1) + O ((1 - z) 2 ) = zA g { 1) = zA« (3.5) 

An extra power of z simply amounts to a shift of N by one unit in the Mellin-space result. 
This clearly does not affect the logarithmic behaviour at N —> oo. This extra z factor 
multiplies also all higher order contributions to the cusp anomalous dimension A(a s ). 

The constants ci have been introduced by requiring 


lim 

TV—>-oo 


Gi(N) - Gi(N) 


= 0. 


(3.6) 


When expanding in powers of a s as in Eq. (3.1), we can effectively get rid of the constants 
by simply defining distributions whose Mellin transform differ by 1/N suppressed terms at 
large N, namely 


1 

i 

i — i 

1- 

+ 

+ 

rin fc 

y/z 

1 

i 

i-H 

1 - z 

1 - z 

1 - z 


(3.7) 


In this way we find 


Gi(A0 = X>? 

n —1 


n n— 1 

E tlMN +1) + E i D lAm 


_k =o 


k =0 


(3.8) 


where the coefficients and j are the same as in Eq. (3.1), and the argument of 
the Mellin transform of the distributions associated with the cusp term are shifted by one. 
Explicit expressions for the Mellin transforms T>i i .(N) are given in Appendix A. Note that 
here, unlike in Ref. [2], we do not shift the D terms; however, the difference is subleading, 
and the impact is negligible. 

It is important to observe that, unlike in the case of Higgs production, the inclusion 
of the second term in the expansion Eq. (3.5), does not lead to the full inclusion of all 
subdominant contributions of the form a™N~ l In 2 ” -1 N. This is because contributions of 
this order to C(N) may arise both from soft emission, but also from interference of powers 
suppressed terms in the leading order cross section. We will use these terms (by turning on 
and off the shift in N) as a way to estimate the uncertainty in our procedure. 
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3.2 Coulomb singularities 


In Sect. 3.1 we have studied the large-IV terms originated by soft gluon emission. In the case 
of heavy quark pair production, however, there is another class of contributions which, order 
by order in perturbation theory, do not vanish in the large-IV limit, altogether unrelated 
to soft emission. It was pointed out many years ago [30] that pair interaction dynamics 
and bound state effects may be relevant in the threshold regime. This multiple exchange 
of “Coulomb gluons” between the two heavy quarks in the final state leads to corrections of 
order ( a s //3) m with f3 = y/1 — p, which are usually referred to as Coulomb singularities. For 
m large enough, such contributions may compete with, and even dominate over, Sudakov 
logarithms in the threshold limit. It turns out, however, that these contributions can be 
resummed to all orders (see for example Ref. [31] and references therein), and the result 
of resummation vanishes in the large-IV limit. Nevertheless, these contributions must be 
taken into account in the absolute threshold limit p —> 1. 

Coulomb terms have been studied in Refs. [29-34], In particular, it was shown in 
Ref. [32] that Coulomb singularities and soft singularities factorize in Mellin space in the 
N —> oo limit and in the singlet-octet basis: 

C(N)= Y, gi(as)MN,a s )e W [G l (N)} + of^-'], (3.9) 

1 = 1,8 ' ’ 


where the factor 

Ji(N,a s ) = 1 + J : (1) (IV)a s + j{ 2) (N)a 2 s + 4 3) {N)a 3 s + 0(aj) (3.10) 

contains Coulomb terms. The coefficients J^\n) and J^ 2 \n) can be obtained by matching 

/o\ 

with exact fixed-order calculations [35], while Jj (N) is unknown. 

It can be shown that the effect of Coulomb terms on the physical cross section is very 
small. In Ref. [29] the impact of the contribution Jj (N) on the hadronic cross section 
was estimated to be about 0.5% of total NNLO correction, and even less at higher orders. 
Similar conclusion can be drawn from inspection of Fig. 1, where we compare, in N space, 
the coefficient function Eq. (3.9) with and without the Coulomb factors Ji(N,a s ), both at 
NLO and NNLO. The effect at N ~ 2, which is the relevant region for Mellin inversion (see 
Section 5.1), is very small. 

A complete resummed expression of Coulomb contributions has been obtained in the 

context of potential non-relativistic QCD (pNRQCD) [36-38]. This calculation can be used 

to test whether the pattern observed at NLO and NNLO continues at higher orders. In 

Fig. 2 we present the same comparison as in Fig. 1 at N 3 LO, with J^(N) extracted from 

pNRQCD calculations [36-38]. The effect of the inclusion of the O^a^) term in Ji(N,a s ) 

(3) 

as predicted by pNRQCD computations is indeed very small. We will include Jj '(IV) as 
computed by pNRQCD methods in our results, but these observations show that our final 
results are not significantly affected by it. An explicit form for all the Coulomb corrections 
is given in Appendix A, Eq. (A. 12) (see also Ref. [39]). 
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Figure 1. The contributions to C(N) from soft gluon emission exp [Gi(iV, a s )] (Sudakov terms) 
and Coulomb terms Ji(N,a s ) (Coulomb terms), at NLO (left) and NNLO (right). The combined 
effect Eq. (3.9) is also shown (Complete soft approximation). 



Figure 2. Same as Fig. 1 but at N 3 LO. The two different Coulomb curves differ by the inclusion 
of the pNRQCD estimate of jj :>> (see text). 


3.3 Final prescription for the soft contribution 

We conclude this section by giving our final prediction for the soft-emission contribution to 
the total cross section to N 3 LO. We define 


CA-softo(-^0 — ^2 gi(a s )Ji(N,a s ) 


1 = 1,8 

( oo 


ex Pi E 


a„ 


, n=l Lfc=0 


E 6 i>‘( w )+E CAw 


n—1 


(D) 


k =0 


CA-softi (-/V) = ^2 9i{oi s )Ji(N,a s ) 


1 = 1.8 


(3.11) 
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exp 


oo 

£ 

< n=1 


a c 


£ b { n y t (N + i )+ y , Ci®* w 


71—1 


(£>) 


,k=0 


k =0 


(3.12) 


where we have made explicit the dependence on a s = a s (m 2 ) and it is understood that the 
exponentials should be expanded in powers of a s up to order a^. We will take the average 
between the two approximations, 

C-A-soft (N) = i [C-A-soft! (N) + CA-softo (N)} (3.13) 

as the central value of our prediction, and the difference 


A soft (lV) = \C A . so{tl (N) - CA-softo (N) I 


(3.14) 


as an estimate of the uncertainty of our procedure. Hence, our soft approximation to the 
coefficient function will be finally given by 

A Soft 

C(N) = CA-soft ± 2 • (3.15) 

We also consider the result which corresponds to the standard resummation, as given 
to NLL in Ref. [33], and extended to NNLL in Ref. [12] (and in the associate public 
code top++), which is expressed as a function of positive powers of In IV, which we will 
call JV-soft, following Ref. [2], This is given by 


Cjv-soft = ^2 9i{®s)Ji(N,a s ) 
1 = 1,8 


exp 


£ 

< n =1 




^ b21v l ° g (N) + 


k =0 


J n,k k 


71—1 

£ 

k =0 


b (D) plog 

°n,k,l U k 


(N) 


(3.16) 


where the functions P[° s (iV) are the large N limit of Vf c (N) written in terms of positive 
powers of In IV, excluding constants (for an explicit definition, see Ref. [2]), i.e. such that 


V k (N)=V^(N) + d k + O^Y 

thereby leading to 


9i(a s ) 


gi{a s ) exp 


S a s ( S }j( r2l d k + J2 b( nh d k 


71—1 


(D) 
k, T 


n= 1 


\k =0 


k =0 


(3.17) 


(3.18) 


4 Small-lV contributions 


In order to extract the leading small -N singularity of the partonic cross section, we will use 
the so-called high-energy or k T factorization technique, first described in Ref. [40] for the 
total cross section, and more recently extended to rapidity distributions [41]. We follow the 
resummation procedure developed by Altarelli, Ball and Forte (ABF) [42]. For a detailed 
derivation of resummation, which affects both coefficient functions and evolution of the 
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parton densities, we refer the reader to the original literature (e.g. [40, 42-44]); here we 

only summarize the main results which are relevant for our discussion. 

In the k T factorization formalism, small-A singularities are obtained by computing 

the leading-order partonic cross section for the relevant process, with off-shell incoming 

2 

gluons. We therefore define an off-shell partonic cross section ^& 0 ff- s heu(p, £i, £ 2)1 which 
is a function of the scaling variable of the process p, and of the transverse momenta of 
the two off-shell gluons: £* = k^/m 2 . Resummed results can be obtained through the 
determination of the so-called impact factor 


h(N,M l ,M 2 ,m 2 ,a s ) = M 1 M 2 R(M 1 )R(M 2 ) 


x / dp p 

Jo 


N—2 




m 
Pv 
— 1 


2 \ M1+M2 


10 


00 o 2 

d£ 1 Zi^^—^Vo&shewiP, £l>6)- 
m z 


(4.1) 


The prefactor R(M\ )R(M 2 ) accounts for factorization scheme dependence [43]; in the MS 
scheme 

(4.2) 


R(M) = 1 + -C 3 M 3 + 0(M 4 ). 

o 

The impact factor Eq. (4.1) for the production of a heavy quark pair was computed in 
Ref. [45]. Here, we are interested in its expansions in powers of M\,M 2 , in the vicinity of 
N = 1: 


h(N, Mi,M 2 , m 2 , a s ) = +0(N-1). 


(4.3) 


11,12 


The coefficients h t] i 2 relevant at N 3 LO can be easily obtained by performing the expansions 
of the previous formula and they are given in the Appendix B, Eq. (B.l). Following Refs. [42, 
46], the resummation is performed by identifying the Mellin variable M* with the resummed 
DGLAP anomalous dimension (together with sub-leading running-coupling effects): 


= 


7 r ; 


where the right hand side is recursively defined by 


+ k+1 
7res 

= 7i(i + *%) 

k 

7res 

L J 

V 7res J 

L J 


[^+ ] = -A 

5 L ^resj /r< 


with 


y + = -/3na 2 -^-7 + 

/ res s dot, res 

By expanding the anomalous dimension to fixed perturbative order 
7 + = a s 7^ + a 2 7 W + a 3 7 ® + O(o.g), 


(4.4) 


(4.5) 


(4.6) 


(4.7) 


we have all the ingredients to construct our A^ —> 1 approximation, which is simply found 
by substituting the expansion (4.7) in Eq. (4.5), and then in Eq. (4.3). The explicit form 
of the anomalous dimensions in the small-lV limit to the order relevant in our discussion 
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is given in Appendix B, Eq. (B.2), (B.3) and (B.4). Now we define a resummed coefficient 
function by factoring out the leading order contribution, evaluated in the high-energy limit, 
i.e. N = 1 


Cabf(A0 — 


at 


m 2 <7 LO (m 2 ,1) 


E h < 

i 21,22=0 


7 r ; 


- 1 


= 2hi,o7 (0) + ( 2/l 2 ,o + * 1 , 1 ) 7 (0)2 - 2h 2 ,o/3o7 (0) + 2/ii,o7 (1) 

(h- 3,0 + h 2 ,i) 27 (0)3 - (3h 3i0 + h,2,i)2/3o7 (0)2 + 4/i3,o/3o7 (0) 


+ O-x 


+ ( 2 h 2 ,o + ^i,i) 2 7^7^ - 4h 2)0 /3o7^ + 2/ii,o7^ + 0(07), ( 4 - 8 ) 


where, in the second equality, we use the fact that 




cf LO (m , 1) = —| h 0 .o, 




(4.9) 


and we dehne 


''* 1,22 


^ 21,22 

ho,o 


(4.10) 


We have omitted the dependence of the coefficients on p 2 for simplicity. 

Since we are going to combine the small -N approximation with the large -N approxima¬ 
tion to obtain an estimate of the full coefficient function, we require that the two limiting 
behaviors do not interfere with each other. In particular, we require that the high-energy 
contribution vanishes when N —> 00 . Manifestly, Eq. (4.8) does not fulfil this requirement, 
because of the presence of a constant contribution in 7 I 0 ) (see Eq. (B.2) in Appendix B), 
which propagates in C'abf(A^) to all orders in a s . Therefore, following Ref. [2], we replace 
Eq. (4.8) with a modified version of the small -N approximation to the partonic cross section, 
which has the same leading singularity in N = 1 but vanishes as N —> 00 . The modified 
version is given by 


C A BFsub(N) = C abf (N ) - 2C abf (N + 1) + C abf (N + 2). (4.11) 


The subtraction simply introduces subleading singularities in N = 0 and N = — 1, but now 


lim C A BF-sub(A0 — 0 (4-12) 

N^-oo 

to all orders in a s . As pointed in Ref. [2], this choice for the subtraction is a compromise 
between the contrasting goals of not changing the small-A r singularities structure and of 
damping strongly enough as N increases. In momentum space the subtraction of Eq. (4.11) 
corresponds to damping the p —> 1 behaviour of the partonic cross section through a 
multiplicative factor (1 — p) 2 . 

One additional modification is needed to dehne our prediction for the cross section in 
the high-energy regime. We note that the anomalous dimensions vanishes at IV = 2 due to 
momentum conservation. This implies that the small -N approximation Eq. (4.8) vanishes 
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in A = 2. This value of A marks the transition between the small-A approximation (not 
accurate if A > 2) and the large-A approximation (not accurate when A < 2). 

However, Eq. (4.8), differs from 0 in A = 2 because the small-A limit of the anomalous 
dimension 7 + ) Eq. (4.7), contains only the leading and the next-to-leading singularities in 
A = 1, and not a full fixed order expression. Momentum conservation can be enforced [42] 
directly in Eq. (4.8), by adding to C ABF (N) a function / mom oc 1/A. Following Ref. [2], we 
construct our small-A momentum-conserving approximation as 

0,e(iV) = C ABF . aub (N) - N ( N k ™^ S \ 2 y (4.13) 

where the constant i m0 m(o s ) is fixed by requiring that Ch. e (2) = 0 to all orders in the 
perturbative expansion. The enforcement of momentum conservation in our small-A ap¬ 
proximation allows us to estimate the contribution of subleading poles in A = 0,-1, —2,..., 
which are not controlled by the high-energy approximation. Indeed, the full partonic cross 
section does not vanish at A = 2: only the contribution to it driven by hard radiation 
from external legs does. Hence we can estimate subleading small-TV effects by allowing the 
partonic cross section to deviate slightly from zero at A = 2 due to contributions from 
subleading poles. We thus fix k mom to be 

^mom(o : s) — C , ABF-sub(2) i ^^mom^s)) (4-14) 

and we take the variation A/c mom («s) to reach as its maximum value 15% of the size of 
soft contribution CA-soft(A), Eq- (3.13) at A = 2. This is a somewhat more conservative 
estimate of the uncertainty in comparison to the corresponding one adopted in Ref. [2], 
due to the fact that we expect the contribution from subleading poles to be potentially 
somewhat larger in this case, based on the behaviour of known orders, and possibly also 
due to the issues discussed in Sect. 2.2. Namely, we choose in Eq. (4.14), 

A k mom (a s ) = 0.15 x CA-soft^). (4.15) 

In conclusion, this means that the small A contribution, rather than being completely 
switched off at A = 2, is small at this point, and gets switched off somewhere in its vicinity. 

5 Approximate cross section up to N 3 LO 

We are now ready to present our results for top pair production cross section at N 3 LO, by 
using Eq. (2.45) to combine the large-A terms Eq. (3.15) and the small-A terms Eq. (4.13). 
We first recall how A-space parton-level results can be related to physical hadron-level 
results using the saddle point methods, and then we present the parton and hadron level 
results in turn. 

5.1 Saddle point approximation of Mellin inversion integrals 

The physical hadronic cross section can be related to the underlying partonic cross section 
by viewing the former as the inverse Mellin transform of its factorized expression Eq. (2.6) 

I /*C - t - 200 ~1 ( t C~\-%00 

(j(m 2 ,p h ) = — dNp- N+1 Jf(N)a(m 2 ,N) = — dN exp[/(A)], (5.1) 

Z7TZ Jq —200 Z7TX Jc—ioo 
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t-tbar - Saddle point N 0 as a function of Vs - m t = 173.37 GeV 



Vs [TeV] 


Figure 3. Position of the saddle point as a function of a/s for fixed top quark mass. 


where c is to the right of all singularities of the integrand. It can be shown on general 
grounds that f(N) has a unique minimum at N = Nq on the real N axis, which allows us 
to estimate the integral by the saddle-point technique. 

The position No of the saddle typically depends very weakly on the partonic cross 
section, and is mainly determined by the parton luminosity; the saddle-point approximation, 
with the inclusion of quadratic fluctuations around the saddle, turns out to be generally 
quite accurate [47]. It is then possible to infer properties of the hadronic cross section from 
the behaviour of the partonic cross section at the value of N which corresponds to the 
saddle point for given hadronic kinematics. 

The position of the saddle point Nq as a function of the collider energy ,/s, for m = 
mt = 173.37 GeV, computed using NNPDF 3.0 [48] parton distributions, is shown in Fig. 3. 
The value of the saddle turns out to be around N = 2.5 at LHC energies. The value of Nq 
is a very slowly decreasing function of the total energy. 

5.2 Parton-level results 

We can now combine the results of Sections 3 and 4 to obtain an estimate of the coefficients 
of the perturbative expansion for the total production cross section of heavy quark pairs. 
We first test our procedure against exact fixed-order calculations, which are available at 
NLO [6, 8] and NNLO [12]. At NLO we specifically compare to the fit to the analytic result 
of Ref. [8] presented in Ref. [49] . In Fig. 4 we show our approximation to the NLO and 
NNLO in Mellin space (labeled approx), compared to the exact results of Ref. [12, 49]. The 
large -N contribution, Ca~so&,(N) Eq. (3.13) (labeled as soft), and the small -N contribution, 
C'h.e(AI) Eq. (4.13) (labeled as high-energy), are also shown. The renormalization and 
factorization scales are taken to be equal to the heavy quark mass m. The uncertainty on 
our prediction (red band) is obtained as the envelope of the uncertainty Eq. (3.15) on the 
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Figure 4. Comparison between exact results and our approximation at NLO (left) and NNLO 
(right) for the Mellin-space coefficient function C(N)\ the large -N contribution (A-soft) Eq.(3.15), 
the small -N contribution (high-energy) Eq.(4.13), and the combined approximation of Eq. (2.45) 
(approx) are shown. The bottom plots show the ratio of the approximate to the exact result. 



Figure 5. Same as Fig. 4, but at N 3 LO. 


soft terms (green band), and the uncertainty Eqs. (4.14), (4.15) on the high-energy terms 
(blue band). 

The agreement is excellent at NLO in the whole range displayed in Fig. 4. At NNLO 
there is a slight discrepancy in the region between 1 < IV < 1.6. This is due to the fact 
we are only including the LL contribution in this region, while it was noted in Ref. [8, 12] 
that NLL contributions close to N = 1 for top pair production are sizable. This region of 
N, however, would only be relevant for very high energy colliders (y/s > 100 TeV). 

At the N values which are relevant for LHC energies, the agreement between the 
approximate and exact results is excellent; for lower collider energy, as the saddle point 
moves towards larger values, the uncertainty on the approximate prediction is smaller. 

We now turn to the N 3 LO contribution. We recall that the constants g^\ together 
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(3) (3) 

with the difference between coefficients Dg ' and D\ , are not known. They will be set to 
zero for the time being; the impact of this missing information on our prediction for the 
physical cross section will be discussed in Section. 5.4 below. The function Capprox(Af) is 
plotted in Fig. 5, together with the soft and the high-energy contributions that go into it 
according to Eq. (2.45). 


5.3 Hadron-level results 


We come now to the hadronic cross sections. It is convenient to define the gluon-channel 
A'-factors 


Kgg(m 2 ,Ph) 


(Jgg(m 2 ,p h ) 

a^(m 2 ,p h ) 


1 + a s K$ + c? a K$ + a 3 a K$ + 0(a 


99 


99 


(5.2) 


where a s = a s (m 2 ), (J gg {m 2 , ph) is the contribution to the hadron-level cross section 
Eq. (2.1) from the gluon-gluon subprocess, and (rri 2 , ph) the corresponding leading-order 
approximation. All results will be obtained using the partonic cross sections of Sect. 5.2, 
with factorization and renormalization scale p, R = p F = m, and NNPDF 3.0 NNLO parton 
distribution functions [48], with a s (M ]|) = 0.118 and rif = 5. Scale uncertainties on our 
final results will be discussed in Sect. 5.4 below. 



Vs[TeV] 


Figure 6. Comparison between the exact result and our approximation for the NLO contribution 
to the A'-factor Eq. (5.2) from the gluon channel, plotted as a function of the collider energy yfs. 
The result obtained expanding the standard NNLL resummation (A-soft) is also shown. 


The NLO and NNLO A'-factors at a pp collider are shown in Figs. 6 and 7, respectively, 
as functions of the center-of-mass energy y/s, and compared to the exact results of Ref. [12, 
49]. We also show the values obtained by expanding out to the given order the standard 
resummed result of Refs. [12, 33], i.e. using the IV-soft approximation Eq. (3.16). The 
main result of this work, namely the N’^LO contribution to the A'-factor in the gluon-gluon 
channel as a function of the collider energy is shown in Fig. 8. Numerical results for the 
A'-factors at LHC energies are collected in Tab. 1. 
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Figure 7. Same as Fig. 6 at NNLO. 
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Figure 8. Same as Fig. 6 at N 3 LO, where the exact result is not known. 


We note that our approximation reproduces the exact result within the estimated un¬ 
certainty both at NLO and at NNLO, in the whole energy range displayed in the plots. In 
comparison to the case of Higgs production, the uncertainty is larger, both because thresh¬ 
old resummation is known to a lower accuracy, and also because, as already mentioned, 
we have less control over ln N N terms. Our result is seen to differ by varying amounts at 
each order from the IV-soft one, simply obtained by expanding out the resummed result. 
At NLO and NNLO the origin of this difference is twofold: first, our approximation also 
includes the high-energy contribution and second, the functional form of the contribution 
obtained by expanding out the iV-soft result differs from that adopted in our approxima¬ 
tion, as explained in Sect. 3.3. At N 3 LO two further differences are, first, that our result 
includes the single logarithmic term, even if its coefficient is only partially known, which 
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a s {m 2 )K 

a 2 s (m 2 )Kjf 

°$( m2 )Kgf 

;> 

exact 

0.599 

0.237 


£ 

approx 

0.623 ±0.112 

0.250 ± 0.076 

0.098 ± 0.038 

L- 

IV-soft 

0.672 

0.260 

0.056 

> 

exact 

0.587 

0.227 


EH 

approx 

0.609 ±0.115 

0.239 ± 0.076 

0.094 ±0.037 

00 

IV-soft 

0.658 

0.250 

0.052 

;> 

exact 

0.555 

0.199 


Eh 

approx 

0.569 ±0.125 

0.209 ± 0.077 

0.083 ± 0.036 

CO 

rH 

IV-soft 

0.619 

0.221 

0.040 


exact 

0.552 

0.196 


Eh 

approx 

0.565 ±0.126 

0.205 ±0.077 

0.082 ± 0.036 

r—i 

IV-soft 

0.614 

0.217 

0.038 


Table 1. The NLO, NNLO, N 3 LO contributions the gluon channel K -factor at the LHC (7 TeV, 
8 TeV, 13 TeV, 14 TeV). 

is absent in a NNLL resummation. However, the numerical impact of this contribution is 
small. Second, that the function of a s which multiplies the resummed exponent is given by 
gi(a s ) in the IV-soft resummed result, and by gi(a s ) when constructing our approximation, 
the relation between the two being given by Eq. (3.18). We will specifically discuss the 
impact of this choice in the next Section (see in particular Tab. 2). 

5.4 N 3 LO top pair production cross section at LHC 

We collect now final results for top pair production at the LHC energies and its uncertainty. 
Our prediction for the total cross section with = /% = m is 

LHC7: u^ pp L r ° x (p h , m 2 ) = (l77.43 ± 2.99 + 0.10 g {3) - 0.10 (d^ - ) pb, (5.3a) 

LHC8: UaSro X (ph-/rn 2 ) = (253.98 ± 4.35 + 0.14 g (3) - 0.14 (d® - D {3) ^ pb, (5.3b) 

LHC13: aJ p L r ° (p h , m 2 ) = (835.61 ± 14.78 + 0.51 g (3) - 0.46 (l>^ 3) - D®)) pb, (5.3c) 

LHC14: aJ p L r ° (p h , m 2 ) = (988.57 ± 17.55 + 0.61 g^ - 0.54 (d (3) - D®)) pb, (5.3d) 

where 

9 {i) =9?+9$ ) - (5.4) 

The constant g ^ is not known. The difference between the coefficients and 
parametrizes the missing information about the single log at N 3 LO, as discussed in Sec¬ 
tion 3.1. Eq. (5.3) is obtained using exact expressions for the NLO [8] and NNLO [9-12], 
combined with our approximation for the N 3 LO in the gluon channel. 

We now consider various sources of uncertainty. First, we consider the uncertainty 
related to missing coefficients. In Tab. 2, we list the values of the coefficients g^\ and 
of D g — D\ for the two known orders, as well as the coefficients defined using the 
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Df - Df g® giW giW - ijW 

a s -0.955 1.188 3.832 2.644 

-1.782 0.535 8.512 7.977 

a 3 ? ? ? 12.348 


Table 2. The coefficients D^' — D^\ g^\ gW and the difference between the latter two at NLO 
and NNLO. 


LHC8 

°qq(P b ) 

Vqgipb) 

a 99 (P b ) 

LO 

29.94 

0 

116.33 

NLO 

4.17 

-0.45 

66.66 

NNLO 

2.17 

-2.30 

26.36 


Table 3. Comparison of contributions to top pair production cross section from different channel 
at LHC 8 TeV at LO, NLO, NNLO. 


analogue of Eq. (5.4) but for the expansion coefficients of gi(a s ) Eq. (3.18) which is used 
in the standard resummed results of Refs. [12, 33], and the difference between the two, 
which is known to a higher perturbative order. Based on these values, we note that the 
perturbative behaviour of the coefficients g® appears to be rather more stable than that of 
the coefficients g^\ and we reasonably expect the unknown coefficients, g ® and — D^\ 

to be both of 0(1). We conservatively estimate the uncertainty related to these unknown 
coefficients as 

5 {3) = 0 ± 5, = 0 ± 10, (5.5) 

which we add in quadrature to the approximation uncertainty of Eq. (5.3). We also note 
in passing that the iV-soft approximation assumes g 0) = 0, while our estimate corresponds 
to a value of order ten for this constant, thereby partly explaining the difference observed 
in Fig. 8 between our result and the IV-soft approximation, as discussed in the previous 
Section. 

A second source of uncertainty is due to the fact that we only predict the contribution 
of the gluon channel. In Tab. 3, we show the contributions to the top-pair production 
cross section coming from different partonic. channels at LO, NLO, NNLO, for the specific 
case of LHC at 8 TeV. We note that gluon fusion contribution is always dominant, and 
largely dominant at high perturbative orders. This remains true for other LHC energies 
(while at the Tevatron the qq component becomes dominant). Therefore, we expect that 
our estimate for the N 3 LO contribution based on the gg channel only is only mildly affected 
by the lack of knowledge of the other channels. However, in order to take into account the 
uncertainty due to the missing partonic channels, we consider the dependence of our result 
on the factorization scale. Indeed, all the scale dependent terms at N 3 LO are available (for 
all channels), but the inclusion of all of them would consistently make the residual scale 
dependence of relative order af. On the other hand, if we include /i F dependent terms 
only in the gg channel, the residual scale dependence is of order ct 3 , since it misses the 
compensation between channels. The difference between these two ways of varying the 
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Figure 9. Factorization (left) and renormalization (right) scale dependence of the top production 
cross section at various perturbative orders. At N 3 LO the factorization scale dependence is shown 
both including the contribution from all channels (NNNLO) and from the gluon channel only 
(NNNLOgg). 

scale can be thus taken as an estimate of the size of the contribution from the missing 
quark channels. 

Finally, uncertainties related to missing higher-order terms can be estimated by varying 
the renormalization and factorization scales in the usual way. The dependence of the cross 
section at LO, NLO, NNLO and approximate N 3 LO on the factorization and renormaliza¬ 
tion scales, is shown in a wide range in Fig. 9; the NNLO+NNLL result of Refs [12, 33] 
is also shown. The factorization scale variation at N 3 LO with fi F is shown both retaining 
contributions from all channels, and from the gg channel only. The factorization scale de¬ 
pendence is rather mild already at NNLO, and even milder at approximate N 3 LO when all 
channels are included. When only the gg channel is included, a stronger scale dependence 
is observed, particularly at the extremes of the range, where sizable contributions from 
the missing channels are generated. The renormalization scale dependence is somewhat 
stronger than the factorization scale dependence, but at N 3 LO it has flattened out almost 
completely, thereby indicating a good perturbative convergence. It is interesting to observe 
that the NNLO+NNLL result has a milder scale dependence than the NNLO, but still a 
stronger scale dependence than the approximate N 3 LO. This is what one may expect based 
on the observation that effectively, because the process is far from threshold, the resurn- 
mation is providing some approximation mostly to N 3 LO, which is however less complete 
than the approximation which is constructed here. 

Our final results, with full uncertainty, are thus 

LHC7: cr^ p 3 j[(o X = 177.43pb ± 1.79%(approx) ± 0.97%(channels)^ 3 '^[(scales) (5.6a) 

LHC8: ^approx = 253.98 pb ± 1.82%(approx) ± 0 . 96 %(channels))( 283 %( sca l es ) (5.6b) 
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LHC13: = 835.61 pb ± 1.88%(approx) ± 0.96%(channels)^2'65%( sca l es ) (5.6c) 

LHC14: CTaj-^ox = 988.57pb ± 1.88%(approx) ± 0.97%(channels)^g2%(scales) (5.6d) 

where the “channels” uncertainty has been computed as (± half) the difference between the 
H F scale variation evaluated with only the gg channel or with all the channels (NNNL0 33 
and NNNLO curves in Fig. 9), in the range m/2 < // F < 2m with /% = m. The “scales” 
uncertainty is instead obtained through a canonical seven-point variation, namely m/2 < 
/%,/r F < with 1/2 < /%/// F < 2, computed with all channels. We observe that the 
approximation uncertainties, though conservatively estimated, are rather smaller than the 
scale uncertainty, and in fact adding in quadrature scale and approximation uncertainties 
we end up with an overall theoretical uncertainty on our N 3 LO result of 3.5%, not much 
larger than the scale uncertainty itself. This uncertainties can be compared to the PDF 
uncertainty, which at the LHC yfs = 13 TeV (with NNPDF 3.0 PDFs) is of order 2%. 
Additional uncertainties come from the values of a s and m;: see Ref. [50] for a more 
detailed discussion. 

We observe that the uncertainty due to scale variations at NNLO is about 5% at 
the collider energies we are considering, which is larger than the overall uncertainty on 
our N 3 LO estimate even accounting for approximation uncertainties. The inclusion of our 
approximate N 3 LO contribution appears thus to be advantageous, and it leads to a decrease 
in theoretical uncertainty which is of the size one would expect when going ftom NNLO to 
N 3 LO. 

We may compare our approximate N 3 LO result to that which would be obtained us¬ 
ing the IV-soft result shown in Fig. 8. The latter leads to a N 3 LO contribution which 
corresponds to an increase of 3%, 2.8%, 2.3%, 2.4% in comparison to the NNLO at LHC 
y/s = 7, 8,13,14 TeV, respectively. This is rather lower than our approximate result, which 
corresponds to an increase of 4.3%, 4.5%, 4.2%, 4.3% respectively. The reasons for this have 
been discussed in the previous section. As discussed above, and as it is apparent from 
Fig. 9, the scale uncertainty on this result is somewhat smaller than that on the NNLO 
result, but larger than that on our result. 

Finally, we note that an approximate N 3 LO result is also presentes in Ref. [51]. This 
result is obtained by considering logarithmic terms enhanced at partonic threshold in the 
differential cross section, inclusive in one of the two heavy quarks produced. The result is 
then integrated to obtain the inclusive cross section. The enhancements over the NNLO 
results were found to be 4%, 3.6%, 2.7%, 2.6% at y/s = 7,8,13,14 TeV respectively. We 
leave a more detailed comparison to the results of Ref. [51], as well as to the resummed 
result obtained in SCET (e.g. Refs. [32, 35] and Refs. [16-18]), for future work. 

6 Conclusions 

We have constructed an approximate expression for the N 3 LO contribution to the produc¬ 
tion cross section of a heavy-quark pair at hadron-hadron colliders. We have focused on 
the gluon-gluon initiated subprocess, which gives the largest contribution at the LHC. 
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We have obtained our result by extending the method developed by some of us for 
the case of Higgs production in gluon-gluon fusion [2], and based on reconstructing the 
Mellin-space partonic cross section from its known singularities. Heavy-quark production 
is a process with a rather more complicated kinematic and color structure compared to 
inclusive Higgs production. Furthermore, fixed-order coefficient function for this process 
are known to have a rather non-trivial singularity structure in physical space, at first sight 
unrelated to the threshold or high-energy limits. Therefore, application to this case of the 
the technique suggested in Ref. [2] and applied there to Higgs production in gluon fusion 
provides a rather stringent test of this methodology. In this study we have shown that the 
method of Ref. [2] provides excellent approximations to known results up to NNLO, thereby 
validating he methodology. 

Having established the reliability of the methodology even in this more subtle case, 
we have used it to produce a N 3 LO approximate partonic cross section for heavy quark 
production. We have then focused on the tt cross section, which is of great interest at 
the LHC. We have found that the approximate N 3 LO correction amounts to an increase 
in comparison to the NNLO prediction of 4.3%, 4.5%, 4.2%, 4.3%, for pp collisions at y/s = 
7, 8,13,14 TeV, respectively. Inclusion of this correction reduces the scale uncertainty to 
3%, with a combined uncertainty on the approximation itself of comparable size or smaller. 
Our final overall conservatively estimated uncertainty is thus somewhat smaller than the 
uncertainty on the exactly known NNLO result, and inclusion of our approximate result 
appears to be advantageous. 
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A Coefficients in the large-A contribution 

The coefficients l) A j, (/i 2 ) and j(^ 2 ) defined in Eq. (3.1) for n < 3 are given by 
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while Dg ' is unknown. 


(A.4e) 

(A.4f) 


To complete our resummed formula for the soft emission, we need explicit expressions 
for the functions D k (N). Mellin transforms of the distributions 


A (z) = 


where the plus distribution is defined by 
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/ dz[f{z)\ + g(z)= dz g (z) [f (z) — f ( 1 )] 
/ 0 Jo 


(A. 6 ) 


The full calculation is presented for example in Ref. [2, 52], Here we give the final result: 
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T(JV) is the usual Euler function, and the superscript in round brackets in T (IV, £) denotes 
differentiation with respect to £. For the first three values of k, relevant for our N 3 LO 
approximation, we find 


A (N) = -L, 
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where ipk(N) = (IV), (N) = r 7 (lV)/r(lV) are the usual polygamma functions and 

7 E = —V’o(l) is the Euler-Mascheroni constant. 

The Coulomb functions Ji(N,a s ) are computed by taking a Mellin transform of the 
resummed momentum-space results obtained in the context of pNRQCD [32, 34]. For 
more details about the procedure of this Mellin transformation, see Ref. [39]. The explicit 
expression for the Coulomb functions, up to N 3 LO, is the following: 


Ji(N, a s ) = 1 + Jj(N)a s + Jf\N)a 2 s + Jf\N)a 3 s + 0{aj) 
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with B (a, b ) = the Beta function. 

We now turn to the functions gi(a s ). They must be computed by matching with exact 
results, which are available up to NNLO. Moreover, as pointed out in Ref. [12], at NNLO 
we can only infer from the exact result the sum g^ = m 2 + g g , but not the individual 
contributions g l . The latter have been estimated by the same procedure adopted in 
Ref. [12], giving 


gi(a s ) = 1 + g[ l) a s + + O (a 3 ) 


(A.13a) 





0.216 


(A. 13b) 
(A.13c) 





+ C A 



0.319. 


(A.13d) 
(A.13e) 


It can be checked that the uncertainty due to this additional guess is completely negligible 
in comparison to our error bands. For simplicity we do not give the explicit factorization 
scale dependence of g(a s ). 


B Coefficients in the small-lV contribution 


The only ingredients which are needed for the high-energy contribution are the coefficients 
h t] of the expansion of the impact factor, up to N 3 LO: 
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The coefficients hi lt i 2 are simply obtained dividing each hi lt i 2 by /io,o- The factorization 
scale dependence can be easily restored following the procedure explained for example in 
Ref. [2], 

The leading and next-to-leading singular contributions to the anomalous dimension 7 + 
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